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Abstract 

Within a variational calculation we investigate the role of baryons for the structure 
of dense matter in the Gross-Neveu model. We construct a trial ground state at 
finite baryon density which breaks translational invariance. Its scalar potential 
interpolates between widely spaced kinks and antikinks at low density and the value 
zero at infinite density. Its energy is lower than the one of the standard Fermi gas at 
all densities considered. This suggests that the discrete 75 symmetry of the Gross- 
Neveu model does not get restored in a first order phase transition at finite density, 
at variance with common wisdom. 
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1 Introduction 


Quantumchromodynamics (QCD) at finite temperature and chemical potential raises 
many interesting questions which may be relevant for the physics of the early universe, 
dense stars or relativistic heavy ion collisions. Since lattice calculations are not yet fea¬ 
sible at p 7 ^ 0, the QCD phase diagram is still to a large extent speculative, although 
there is little doubt that the phase structure is quite rich [|TJ. In view of this situation, 
it may be of interest to go back to some exactly soluble held theoretic models and try 
to understand their phase diagram, in particular the structure of dense, baryonic matter 
(for a recent review, see Ref. |2|]). The simplest such model is perhaps the Gross-Neveu 
model |Q, N species of massless fermions with a quartic self-interaction, 


£ = + —g 2 ('ip'ipf + A^iTsVO 


A = 0,1 . 


( 1 ) 


(Flavor indices are suppressed, i.e., etc.) This model can be solved 

analytically in the large N limit and shares some properties with QCD, notably asymptotic 
freedom. It exists in two variants with either discrete (A = 0) or continuous (A = 1) chiral 
symmetry. For the sake of simplicity, we shall refer to the first model as Gross-Neveu 
(GN) model, to the second one as 2-dimensional Nambu-Jona-Lasinio (NJL 2 ) model Q. 
Aside from the spectrum which contains massive, unconfined fermions as well as composite 
mesons and baryons, the phase diagram of these models has also been studied some time 
ago [§ § 0- 0’ 0- 0]- It is quite non-trivial, including a tricritical point separating lines 
of first and second order transitions in striking analogy with recent conjectures about 
real QCD El . Here, we shall be mainly interested in the behavior of these models at 
finite density but zero temperature. There, either the discrete or the continuous chiral 
symmetry which are spontaneously broken in the vacuum are believed to get restored in a 
first order phase transition at the critical density p cr ; t = m/( 7 rv / 2 ), with m the dynamical 
fermion mass in the vacuum. Below this density, matter should consist of droplets of 
chirally restored phase containing extra fermions surrounded by regions of chirally broken 
vacuum, much like in the bag model. Above p cr i t the fermions become massless and free. 

This scenario which dates back to the mid 80’s was recently challenged in the case of the 


NJL 2 model by one of us [12 , and we extend our criticism to the GN model in the present 
work. Our main objection is very simple indeed: The conventional wisdom ignores the 
fact that these models possess baryons, either massless (NJL 2 model) or massive ones 
(GN model) (see || and references therein). In the low density limit, one would expect 
that baryonic matter will form widely separated baryons, its energy density satisfying 


dE_ 

dp 


= Mi 


B 


( 2 ) 


p =o 


This relation is violated by the standard solutions to the GN or NJL 2 models. The 
reason can be traced back to the tacit assumption of unbroken translational invariance 
inside and outside the droplets. In the case of the NJL 2 model, the energetically favored 
structure of baryonic matter was shown to be a chiral crystal |[l2j : Baryons are topological 
excitations of the pion held, analogous to the Skyrrne model in higher dimensions ||l3| , 
and matter is a kink crystal, the low dimensional analogue of a Skyrrne crystal |14|. This 
result can be shown analytically with the help of bosonization techniques; Eq. (J2j) is then 
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indeed satisfied for the true ground state with M B = 0. Strictly speaking, no restoration 
of chiral symmetry takes place as a function of density, although a rapidly fluctuating 
condensate may be considered as equivalent to the chirally symmetric phase for most 
practical purposes. The same argument will apply for the crystalline configuration of the 
GN model constructed in this work. 
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The problem in the GN model is more involved, since the bosonization used in Ref. 
can only be applied to the light “pion”. Therefore, we cannot present a solution to the 
finite density GN model (possessing only the massive “sigma” meson) which would be 
nearly as complete as for the NJL 2 model. We have performed a variational calculation 
in order to find out whether the standard Fermi gas ground state of baryonic matter in 
the GN model is unstable against crystallization. 


We finally remind the reader that quite generally, breaking of continuous symmetries and 
Goldstone bosons are not allowed in 1+1 dimensions ]16], [17]]. As is well known, the large 
N limit suppresses the fluctuations which would otherwise restore the symmetry. This 
makes it possible to study such effects in 1+1 dimensions as well, provided one works only 
in leading order in 1 /N, i.e., takes the limit N —> 00 at face value. This remark applies 
both to the continuous chiral symmetry of models like the NJL 2 or ’t Hooft model [T3| 
where massless, pion-like excitations appear, and to translational invariance in the GN 
model, our main theme here. 


2 Construction of a trial wave function 


In the large N limit, the relativistic Hartree-Fock approximation for the ground state 
becomes exact. It corresponds to a variational calculation where the trial wave functions 
are all possible Slater determinants; the Hartree-Fock method determines the best fermion 
single particle orbits. Since we cannot solve this problem fully at this point, we search for 
the best Slater determinant in a restricted subspace only, generating the single particle 
orbits from a certain one-parameter family of periodic scalar potentials. These scalar 
potentials need not be self-consistent since they are only used to restrict the subspace 
of Slater determinants. By evaluating the expectation value of the exact Hamiltonian H 
in this trial state and minimizing it with respect to this one parameter, we obtain an 
upper bound for the ground state energy. If this upper bound is below the energy of the 
translationally invariant Hartree-Fock solution, we can in this way establish instability of 
the system with respect to breakdown of translational invariance. In the present section, 
we outline the construction of our trial wave function. It is motivated primarily by our 
expectation that matter at low density will tend to form isolated baryons. Another guiding 
principle is the desire to gain as much analytical insight as possible. 


Our starting point is the scalar potential for a single baryon in the Gross-Neveu model. 

( 3 ) 


Using units where m — 1 (m is the physical fermion mass), it is given by [ 18, 19 


S(x) = 1 + ?/(tanh£_ — tanh£ + ) , £± = yx ± -artanhy . 

The mean-field Dirac Hamiltonian with this potential, 

h D = 75t^- + ^°S(x) , 

1 ox 


( 4 ) 
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can be diagonalized analytically. The parameter y G [0,1] is related to the occupation of 
the valence level. If y > \/3/2, S(x) crosses zero at the points x = ±R/2 with 


J_ j ( V ~ i + yw \ 
2y 11 \2y 2 - 1 - yw) 


w = 


V V - 3 • 


(5) 


Furthermore, we note that as a result of spontaneous breakdown of the discrete 75 symme¬ 
try, S(x) and —S(x) are equivalent scalar baryon potentials built on the two degenerate 
vacua. The corresponding single particle solutions of the Dirac equation can be related 
by il>'(x) = vy^ix). This enables us to construct a family of smooth, periodic potentials 
for which the Dirac equation is again exactly solvable: Choose S(x) from Eq. (^) in the 
interval [— R/2,R/2] and extend it to the whole axis via S(x + R) = —S(x). By pasting 
together pieces of the single baryon potential with alternating signs in this manner, we 
obtain a continuous and differentiable scalar potential whose period of 2 R is controlled by 
the parameter y. The shape of the potential changes with the period in a prescribed way, 
see Fig. 1. In the limit y —> 1, the potential corresponds to a succession of equidistant, 
alternating kinks and anti-kinks at a distance growing logarithmically with (1 — y). In 
the opposite limit y —* y/3/2, S may be interpreted as arising from strongly overlapping 
baryons with an almost complete cancellation of their scalar potentials. Only some rapid 
residual oscillation around zero is left which goes over into the value zero of the chirally 
restored phase at the endpoint y = \/3/2. The resulting scalar potential is not only peri¬ 
odic with period 2 R, but also invariant under translations by half a period (R) combined 
with a discrete 75 -transformation. In this sense, it may be viewed as discrete version of 
the “chiral spiral” of Ref. 0 - What we have set up here is evidently a kink-antikink 
crystal, as opposed to the kink crystal characteristic for the ground state of models with 
continuous chiral symmetry. 


The Dirac equation for a single baryon can be solved in closed form [ITR pi]. Since we shall 


be mostly interested in negative energy continuum states in what follows, we recall the 


well-known results for these particular wavefunctions only. In the interval [-R/2, R/2], 
one finds (we use a representation where 7 0 = — 07 , 7 1 = io- 3,75 = 7 ° 7 1 = —cr 2 ), 


i>k{x) 


(i k - l)(i k -y tanh£_) \ ikx 
-E(k){ik -ytanhf+) ) 


E(k ) = y/W+l . 


( 6 ) 


ipk and are degenerate with eigenvalue —E{k). The general solution of the stationary 
Dirac equation to — E(k ) in this interval is then 


i>{x) = Ai> k {x) + Bi/j_ k (x) . (7) 

Furthermore there exists a discrete valence level at energy \/l — y 1 . 

Let us now apply the textbook procedure from solid state physics in order to extend the 
solution to the whole x-axis and determine the spectrum of the Dirac Hamiltonian. In a 
slight modification of Bloch’s theorem, we require the following boundary condition, 

VW 2 ) = e i< 5 i 7 5 ' 0 (—i?/ 2 ) . ( 8 ) 


This guarantees a smooth matching of the wave functions at the borders of different 
intervals. Inserting the ansatz (j7|) and denoting the upper (lower) components of by 
ip k ( Xk), Eq. (H|) yields the generalized eigenvalue problem 

Mv = e i5 Nv (9) 
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with 


( Mm) <e-k(R/i) \ N = ( -»(--R/2) -xM-m) i „ = M \ 

1 x*c-r/ 2) x_*c«/2) i' v m-m) <p-t(-m))' \ B )' 

( 10 ) 

Hence the factors e 1<5 are the eigenvalues of the 2 x 2-matrix N 1 M. Evaluation of this 
matrix using Eq. (p|) shows that it has unit determinant, whereas its trace is given by 

tr (iV -1 M) = 2X k , (11) 


with 


X k 


-m 

k(k 2 + y 2 ) 


(k 2 + 1 


The characteristic equation then reads 


y 2 ) sin (kR) + kw cos (kR) 


( 12 ) 


(e i<5 ) 2 - 2X k e iS + 1 = 0. 


(13) 


For |X(/c)| < 1, the matrix N~ X M is unitary, 5 = arccos X k is real and the corresponding 
level lies inside an allowed band. For |X(/c)| > 1, 5 is imaginary and we set 6 = i£. The 
eigenvalues of 1V _1 M are now e^ with wave functions blowing up in one direction — the 
corresponding state is forbidden and lies inside a band gap. The diagonalization of Eq. 
(|9j) also determines the ratio A/B of the coefficients introduced in Eq. ( 0 ); one finds 

I = rq? ( Yk ± +? - (i - WJeW) . (i4) 

with 

Y k = kwsm(kR) — (k 2 + 1 — y 2 ) cos (kR) . (15) 


The two solutions correspond to the two degenerate states with energy —E{k). This 
procedure yields both the spectrum of the Dirac Hamiltonian and the single particle wave 
functions up to the normalization for which we shall use a box normalization later on 
(cf. Sect. 3). Here, let us briefly illustrate the resulting band structure of the potential 
used to generate our trial wave function. In Fig. 2, we show an example of the single 
particle spectrum, comparing it to the one corresponding to the interacting vacuum (or, 
equivalently, free massive fermions with m = 1). It exhibits the expected band structure. 
Since all eigenvalues still occur in pairs ±E(k), the Dirac sea inherits this band structure, 
a novel feature as compared to a non-relativistic solid. Fig. 3 displays the lowest few 
bands as function of the parameter y. The left-hand part of this figure (corresponding to 
small lattice spacing 2 R) can be understood perturbatively, since the scalar potential is 
weak (cf. Fig. 1). Using degenerate perturbation theory, one easily finds the following 
position and width of the gaps, expressed in terms of q — y — y/3/2, 


~ (2n + lK / _ _y _ N 
gap ~ { 2 ^ 3 , 


r « 

1 gap ~ 


64y/3t] 
7r 3 (2n + l) 3 


n = 0,1, 2... 


(16) 


On the right-hand side of Fig. 3 (large lattice spacing), all gaps move into the mass 
gap and we recover the spectrum of the Callan-Coleman-Gross-Zee kink [[TE|], including 
the valence state at zero energy. The lowest allowed band obviously develops out of the 
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valence level of the single baryon (which moves to zero energy as y —> 1), see the dashed 
curve in Fig. 3. Finally, in Fig. 4, we illustrate the band structure together with the 
corresponding scalar potential, for three values of y. Case a) features well separated 
baryons, weak tunneling and hence a narrow valence band. In case b), the band width 
is comparable to the depth of scalar potential. In case c), the valence band has become 
much wider than the amplitude of the scalar potential. Fermions can move freely through 
the whole crystal and baryons have lost their identity. By varying the parameter y, we 
let the system decide which option it prefers at a given density. 


3 Computing the regularized ground state energy 


In our variational calculation of the ground state, we proceed as follows: We take the 
periodic scalar potential from the preceding section and fill the single particle states 
until the desired average baryon density is reached. Next, we evaluate the expectation 
value of the exact Hamiltonian in this trial state and minimize with respect to y, our 
only variational parameter. At the minimum, there will be a unique relationship between 
lattice spacing 2 R (or, equivalently, y) and baryon density. In the case of the single baryon, 
y was linked to the occupation fraction v = n/N of the valence level via y = sin(7rn/2). 

Evaluation of the ground state energy cannot be done naively for several reasons. First, 
there is a trivial quadratic divergence in the kinetic energy due to the Dirac sea. This will 
be removed by subtracting the ground state energy of a free Fermi gas of corresponding 
density. Then, there are logarithmic divergences both in the kinetic and in the potential 
energy which require the usual renormalization of the theory at zero density. Moreover, 
due to the gap structure of the Dirac sea, a direct calculation in the continuum would 
be technically difficult. We therefore enclose the system in a box of size L (an integer 
multiple of the lattice spacing, L = 2 RK) and require antiperiodic boundary conditions, 
thereby discretizing the spectrum of hp. The size of the box will be made large enough 
so that our results do not depend on the IR-regulator L. The phase factor e 1<s in Eq. (||) 
then has to satisfy 


(—e 2l5 ) A = —1 , 5 = p n R+^, p n = ^-(2n+l) (17) 

with n an integer. Solving the characteristic equation ([[3]) yields 


e 15 = X k ± i\/l - XI 


(18) 


or, equivalently, 


- sin (p n R) = X k . 


(19) 


A solution of Eq. (pffi|) for k then yields the single particle energies ±E(k) corresponding 
to the discrete p n . 


To reach a certain baryon density, one has to fill all negative energy states plus a number 
of positive energy states. “Antimatter” on the other hand is obtained by leaving a number 
of negative energy levels unoccupied. Since the energy density is identical in both cases, 
we have found it more convenient to lower the level of the Dirac sea rather than raise 
it. This is the reason why — as mentioned above — we need to consider only negative 


6 



energy single particle states below the “Fermi surface”. We denote the single particle 
wave functions by where n refers to the energy and % — 1,2 labels the two degenerate 
solutions. 


Taking the expectation value of H in our trial state and invoking the large N limit yields 
the following expression for the ground state energy per flavor, 

E = T + V (20) 


with 


occ 2 

n i=l 
rR/2 


rR/2 

-R/2 

2 


1 d 


^- 77-^75 

Mn,i 1 OX 


occ 2 


V = —Ng z K I da; MY , 

'-R/2 


W) = E E 


^ A f ' 

l n,i 


The normalization factors are 



( 21 ) 

( 22 ) 


(23) 


The kinetic energy has both quadratic and logarithmic divergences. Let us isolate the 
divergences by adding and subtracting the asymptotic form (n —> oo) of the terms in Eq. 
0 ), 

T = T conv + Tdiv , (24) 

with 


T 

-L C ( 


= e \± ™ n dx - - (i - f)}, 

n J ~ R / 2 N n .i 1 OX p n \ R ) ) 


T, 


div 


- 5Z i + ^ ( 1 R 


Pr 


2w 


(25) 


ffere, p n are the momenta of the free fermions introduced in Eq. dUP- One can derive this 
asymptotic form using the known single particle wave functions of Sect. 2. The potential 
energy has a logarithmic divergence which can be exposed as follows. First, split the 
condensate into regular and singular pieces, 


with 


(V’V’) = (^) conv + (^^)div , 


(V^)conv 


occ 2 


EE 

n i =1 


1 

a r n ,i 







(26) 


(27) 


Once again, the explicit form of the single particle wave functions has been used to derive 
Eq. 0, where S is the scalar potential described above. As usual, the gap equation at 
zero density is used to renormalize the coupling constant, 


1 

iv 


2 ^ 1 
L \ E(p n ) ' 


(28) 
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Here, the summation over n is according to the vacuum occupation, as opposed to the 
occupation for baryonic matter which enters the sums in Eqs. ([HJ22). The UV-cutoff A 
has to be the same in both sums. We then insert the decomposition (|26|) into Eq. (|22|). 
Since Ng 2 vanishes like (In A) -1 , we may drop the ('0'0)conv term and have to keep only 
the two remaining terms, 

V — Vi + V 2 


(29) 


with 


V 1 = —Ng 2 K 
V 2 = —Ng 2 K 


rR/2 

l-R/2 

rR/2 

l-R/2 


d(C 2 (/Ip'lp) conv {jP'lp} div ) 

is (<W>)di v - 


(30) 


Using Eqs. (|7]) and ([281) and performing the limit A —> 00 , V\ reduces to 

r R/2 


V 1 = 2K I dx(^) mnv S. 

l-R/2 


V 2 is given by 


U 2 = 


( Y~>OCC 1 
E n Pn 


rR/2 


(y vac _1 ) R J-R/2 

V^ n E(p n )J 


d xS 2 


(31) 


(32) 


The integral can be performed in closed form with the compact result 

r R/2 
l-R/2 

Assuming the cutoff A to be large then yields 

( OCC 2 


rR/2 

/ d x S 2 = (R — 2 w) . 

J-R/2 


/ vac 1 


n Pn \ n E (Pn ) 


OCC -\ \ 

E- 

n Pn, 



(33) 

-!)• 

(34) 


If we add up T div , Eq. ( |25| ) and V 2 , Eq. ([H]), and eliminate the trivial quadratic divergence 
by subtracting the free Fermi gas kinetic energy at the same density, we obtain the UV 
finite result 

2w\ 1 1 


T div + ^2 + E ^Pn — (1 — 


R J 


E 


E(p r , 


V 

n Pn, 


(35) 


After some rearrangement, the total regularized energy can be cast into the form 
Er, 


reg 


= £ (E 2 A- r + 2 p„ - i (l - f) } 

n [ i^l J - R /2 Nn,i 1 OX p n \ R J ) 


r R/2 


occ 2 


1 


5 


+2K / dx S E "T7 '$n,i'lftn,i + T E 


l-R/2 


i=1 Nn,i 


-b E - 

E (Pn) n Pn, 


(36) 


Invoking the fact that the ip n) i are eigenfunctions of the Dirac Hamiltonian (|4|) with 
eigenvalues —E(k n ) as well as the result (|53|), this expression can further be simplified to 


o 7/ ,\ vac 1 

E ” S = ~ 2 T.( E ( k «)-Pn) + ( 1 " -R ) £ E/y 


(37) 












The k n are the single particle momenta to the corresponding discrete p n , obtained via the 
transcendental equation m- 

Remarkably, the total energy has been expressed entirely in terms of the single particle 
energies and known functions of y. That the sum in E reg is UV finite can now be confirmed 
by means of the asymptotic behavior 




Pn ~ 


w 


RPn 


n 


oo 


(38) 


which in turn follows from the transcendental equation (|l9| ). Summarizing, the calculation 
of the variational energy has been reduced to solving a transcendental equation for k n , 
Eq. (|19|) , and performing a convergent sum, Eq. (|37]). The baryon density is adjusted via 
the occupation of the single particle states. There is only one variational parameter y - 
it controls the spatial period of the scalar potential underlying our trial wave function. 


4 Results of the variational calculation 


We have performed calculations along the lines outlined above and evaluated the ground 
state energy of baryonic matter with our variational ansatz, the kink-antikink crystal. 
The size of the box L was increased until the results were stable which required values of 
L ss 300 — 1600 (in units where m = 1) increasing with density. The sum to be performed 
is UV finite so that no cutoff is required. To get precise results, we had to take into 
account single particle orbits with energies up to \E\ « 50 — 100 typically. 

We first summarize our Endings concerning the variational parameter y. The result is 
extremely simple and can be stated as follows: At all densities considered, the best choice 
of y was the one which yielded a fully occupied valence band in baryonic matter (or empty 
in “antimatter”), cf. the shading in Fig. 4. With hindsight, this is rather plausible: 
Around each gap, the levels are pushed apart symmetrically due to the interaction (in 
first order perturbation theory). When summing over occupied states the main effect 
cancels out. Only at the Fermi surface the filling is asymmetric and the system can take 
advantage of the level shifts to lower its energy. This picture was confirmed by studying 
in detail where the bulk of the energy difference came from, ft was indeed the region 
in the vicinity of the Fermi surface. The optimal value of y translates into the following 
relation between baryon density and baryon radius (or lattice spacing), 


Pb - 


1 

2 R ' 


(39) 


By this relation, the different values for y in Fig. 4 correspond to baryon densities a) 
p = 0.12, b) p = 0.3, and c) p = 1.0 . Each potential well may be associated with one 
baryon, although for high densities fermions are not confined to one well, cf. Fig. 4. 

With this optimal choice of y we get a unique answer for the best ground state energy 
as a function of density. It is shown in Fig. 5, together with the energies of the free 
Fermi gas and the energy for the standard “mixed phase” solution. We recall that in 
the latter phase, the extra fermions are put into droplets of chirally restored phase, like 
in the bag model, whereas the vacuum in between the baryons remains in the broken 
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phase. The standard first order phase transition occurs at pf = or p ps 0.225, as 

indicated by the dot in Fig. 5. Our result is always below these other two curves. At small 
densities, it perfectly matches the expected slope, Eq. (^), with M B — 2/tt, the mass of 
a single baryon with fully occupied valence level. This is no surprise since our variational 
ansatz was designed to have this limit correctly built in, but it is a good test of the whole 
calculation, including the regularization and renormalization procedure. In our ansatz, 
we had to assume that the distance between two baryons equals the baryon diameter in 
order to generate a smooth periodic potential by gluing together pieces of the single baryon 
potentials. This is certainly the right thing to do at very low densities. There, the extra 
fermions are concentrated in the kinks and anti-kinks. The interaction between kinks and 
anti-kinks is repulsive, therefore they will tend to spread equally; there is no distinction 
between left and right neighbours. It is more surprising that such an ansatz is also capable 
of lowering the energy at higher densities, even at densities where the standard picture 
predicts a transition to the chirally restored phase. Thereby, this standard first-order 
phase transition in the GN model can be ruled out, because our variational upper bound 
for the ground state energy excludes the chirally symmetric solution. In a sense, the mixed 
phase picture may be regarded as a crude attempt to describe baryons. It misses some 
important binding effects which survive the large N limit [^(J , as manifested through the 
existence of single baryons in the spectrum of the GN model. 

We have pushed the calculation up to higher densities without any indication that the 
curves shown in Fig. 5 will cross each other, see Fig. 6. The energy difference becomes 
tiny and the volume one needs to get the same accuracy becomes huge, so that it gets 
increasingly difficult to get precise numerical results. Our trial state definitely becomes 
worse with increasing density, as can be judged by comparing the chiral condensate tppif) 
computed from our wave functions with the input scalar potential. Hence the true ground 
state must be lower in energy at high densities than suggested by our variational calcula¬ 
tion. 

Finally, in Fig. 7, we illustrate the spatial dependence of the baryon density for the three 
values of y already used in Fig. 4. Unlike in the NJL 2 model where the baryon density 
was constant, here it also reflects the periodicity of the scalar potential. 

Needless to say, it would be of quite some theoretical interest to find the truly self- 
consistent periodic Hartree-Fock potential for the GN model at finite density. This may 
be the prerequisite for revising the phase diagram of the GN model in the whole (T, p) 
plane. 
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Figure 2: Single particle spectrum in a reduced zone scheme as function of the Bloch 
momentum p = ($ ~ §), restricted to the first Brillouin zone • Solid line: 

single particle energies obtained from the transcendental equation S(k) = arccos Xj? for 
y = 0.95; dashed line: single particle energies in the interacting vacuum. 
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Figure 4: The lowest band gaps (straight, solid lines and shaded areas) in comparison to 
the corresponding scalar potential (wiggly line), for a) R = 4.16, b) R = 1.66, and c) 
R = 0.5. The hatched regions are the occupied levels in the antimatter crystal. 
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Figure 5: Energy density per flavor as function of the baryon density pB- Dashed line: 
non-interacting Fermi gas; dash-dotted line: standard mixed phase configuration, with 
the first order phase transition at the dot; solid line: variational upper bound for the 
kink-antikink crystal. 
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Figure 6: Difference of the energy densities per flavor between the crystal and the trivial 
Fermi gas solutions. The negative sign shows that the crystal is energetically favored at 
all densities considered. 
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Figure 7: Spatial dependence of the baryon density (dashed), compared to the scalar 
potential (solid), for mean baryon densities a) p = 0.12, b ) p — 0.3, and c) p = 1.0. For 
p = 1 the oscillation of the baryon density is too weak to be resolved. 
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